Sparse repulsive coupling enhances synchronization in complex networks 
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Through the last years, different strategies to enhance synchronization in complex networks have been 
proposed. In this Letter, we show that the synchronization in a small- world network of attractively coupled 
non-identical neurons is strongly improved by adding a tiny fraction of phase-repulsive couplings. By a 
purely topological analysis that does not depend on the dynamical model, we link the emerging dynamical 
behavior to the structural properties of the sparsely coupled repulsive network. 
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Synchronous oscillations in a large ensemble of oscilla- 
tors are considered as one of the mechanisms in biological 
networks to transmit and code information, especially in 
the brain Q]]. Recent experiments have pointed out the im- 
portant role that the complex structure of connectivity has 
in this collective behavior \2], obtaining the signature of 
an underlying small- world (SW) network by indirect mea- 
sures in neuronal culture samples yB or using functional 
magnetic resonance imaging in humans 0]. Theoretically, 
several strategies have been developed with the aim of find- 
ing the best way to achieve synchronization in complex 
networks |5l@.Q|. These approaches have mainly focused 
on the role that weighted links play in heterogeneous net- 
works 1 8, 9. 11(11 . shortest paths between nodes and cluster- 
ing structure in SW networks 1 1 1 ], or the input degree each 
node receives regardless of the net structure I12ll . 

Most of this research has been devoted to attractively 
coupled dynamical elements. However, it is known that bi- 
ological networks combine different types of connections 
to improve synchronization and transmission performance, 
as in the case of the coexistence of excitatory and inhibitory 
synapses in the brain H 1 311 - Nevertheless, little attention has 
been paid to the effect of repulsive coupling, or to the inter- 
action between different types of coupling. The scarce lit- 
erature addressing synchronization in repulsively coupled 
oscillators considers global or local coupling 1 14l ll5ll . but 
the influence of the network structure is still an open ques- 
tion. 

In addition, almost all the published work on synchro- 
nization in complex networks basically deals with arrays 
of identical units. However, heterogeneity of elements is 
an inherent feature present in natural systems which can be 
especially relevant in the dynamics of biological networks 
llql . In this Letter, we explore the influence of the net- 
work topology on the dynamics of non-identical coupled 
units, when a small fraction of the links is phase-repulsive. 
We show that sparse repulsive links in a SW structure can 
induce a coherent oscillatory state when the equivalent SW 
composed of only attractive connections is not able to syn- 
chronize or even to activate the ensemble. Then, just by 
means of an analysis focused on the connectivity matrix, 
we link the emerging dynamical behavior to the structural 



properties of the sparsely coupled repulsive network. 

We study the dynamics of an ensemble of non-identical 
locally coupled Hodgkin-Huxley (HH) neurons 11711 con- 
sidered as spatially isopotential cells 

CV t = Li — ifiV^Xi) + d^2 LijVj (1) 

3 
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Here, Vi is the voltage across the membrane of neuron % 
of capacitance C, and x = {m, h, h} describes the gat- 
ing of the ion channels. I\ on = gNa^hiiVi — V Na ) + 
QK^iVi— VK")+<7i(Vi— Vi) is the ionic current mainly car- 
ried by Na + and K + ions and other ionic currents through 
voltage dependent channels. These currents are driven by 
the voltage difference with respect to the equilibrium po- 
tentials Vjva, Vk and V/ and the maximal ionic conduc- 
tances g Na , g K and <?;. Functions a x and j3 x are voltage- 
dependent rate constants. Parameter values and functions 
are the standards in the literature Ii accounts 

for any external bias current, which has been chosen as a 
control parameter to introduce heterogeneity in the popu- 
lation by setting Ii uniformly distributed within the inter- 
val I ± AI. The value I = 9 //A/cm 2 is fixed close 
to the point where an inverse Hopf bifurcation occurs, 
lb = 9.5 /uA/cm 2 . This way, for the chosen AI = 0.2, 
about 90% of the neurons stay around the silent state, 
V KSt = —65 mV, while the rest will fire periodically. 

The coupling structure in ©is given by = Ly/fc,-, 
where is the Laplacian matrix 11911 . fc, normalizes the 
connection strength by the number of incoming links to 
node i, and the coefficient d stands for the global coupling 
strength. 

Local coupling. — Initially we consider the effect of a 
regular lattice on an ensemble of N neurons, both for fully 
phase-attractive and phase-repulsive coupling. The Lapla- 
cian matrix for the first case is -L^±i = 1, La = —2 
and L 

Figure \l\ shows the global mean firing rate (MFR) and 
its standard deviation ct M fr as a function of d ranging from 
negative to positive values. The negative sign of d comes 



ij — otherwise. And, for the second one, it is 
1, L„ = 2 and L, , = otherwise. 
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FIG. 1 : Mean firing rate (left) and its standard deviation (right) 
for a TV = 400 regular array as a function of d, each point aver- 
aged over 100 realizations. The negative sign of d comes from 
the fully phase -repulsive connection matrix. 



from the fully phase-repulsive connection matrix. When 
d > is large enough, the system is frequency entrained 
to a phase synchronization state. Equivalently, for a suf- 
ficient d < 0, the system reaches an anti-phase synchro- 
nization state. It can be noted from Fig.^that the entrap- 
ment with negative couplings is achieved for smaller abso- 
lute values of d compared to the case with positive ones. 
This indicates that a phase-repulsive coupling is more ef- 
fective to activate and entrain the whole network. Many 
biological systems exhibit this kind of repulsive coupling 
when their dynamical units are in competition with each 
other. Known examples are the inhibitory coupling present 
in neuronal circuits associated to a synchronized behavior 
in central pattern generators [20] or calcium oscillations in 
epileptic human astrocyte cultures 12 ill . 

Non local random coupling. — Our main interest is to 
explore the influence of a SW-like connection topology in 
the activation and synchronization of the network. From 
the results obtained in the previous section, we know that 
a small positive coupling strength is less efficient than a 
negative coupling to activate and synchronize the whole 
network. Taking this into account, we consider now the 
possibility of non-local couplings both positive and repul- 
sive. The global coupling strength is fixed to d = 0.1, i.e., 
within the unsynchronized regime for local positive cou- 
pling as shown in Fig. The Laplacian matrix L is mod- 
eled now by keeping the regular short-range connections 
positive, L iA ±i = +1, and by randomly adding (rather 
than rewiring) a fraction p of the (N — l)(iV — 2)/2 pos- 
sible long-range links, being negative with probability p n . 

Figure |3 shows space-time plots of the voltage variable 
through the whole array for different probabilities p and 
p n . As expected, in the absence of long-range connections, 
few more than the initial 10% of the neurons is firing for 
the chosen coupling strength d, i.e., the array is not even 
activated as shown in Fig. a). When long-range links are 
included, the first observation is that for any p, a minimum 
fraction of the new added links needs to be repulsive in or- 
der to increase the activity of the network. This becomes 



evident when comparing Fig. 0b) with Figs. 0c)-(e). In 
Fig. 0b) the activity generated by the 10% of initially ac- 
tive neurons is reduced, or even annihilated, when all long- 
range connections are attractive (p n = 0). However, the 
scenario completely changes when, for the same p, some 
of the shortcuts are repulsive (p„ > 0) as in Figs. 0c) -(e) 
where self-sustained electrical activity emerges for nonzero 
p n . In addition, we observe the existence of optimal prob- 
abilities p and p n for which the collective oscillation be- 
comes maximally phase-coherent. This fact can be ob- 
served by comparing Fig. 0c), where p and p n are optimal, 
and Fig.0d) where p n is the same but p is slightly higher. 

To study quantitatively how the dynamics is affected by 
p and p n , we measure the MFR of the network and the 
standard deviation of the global electrical voltage, V(t) = 

Eti Vi(t), obtained as a v = ^WWTWW- 
While the MFR gives us an estimation of how much the 
network is activated, the oy defines how coherent is the 
activity of the entire network. If the network is fully ac- 
tivated, the MFR approaches to a rate of around 70 Hz, 
whereas ay is maximal if this activity is synchronized. 

We have plotted in Fig. both the MFR and the ay as 
a function of the probability p for different values of p„. 
The signature of a network resonance both in p and p n is 
clear from this figure: the frequency entrainment increases 
and phase synchronization is maximally enhanced for the 
optimal values p = p c and p n ~ 0.3. The probability p c 
depends slightly on p n , shifting to higher p as p n increases, 
but remaining very small. The interplay between topology 
and dynamics becomes evident when we observe that the 
critical link probability depends strongly on the size en- 
semble as p c oc \n(N)/N, i.e., coincides with the birth 
of the giant connected component (GCC) when only the 
random part of the network is considered. 

Recently the method of the master stabil- 

ity function 123B has been successfully used to analyze 
whether the network structure has some bearing on the dy- 
namics evolving on it. However, this approach requires the 
dynamical units to be identical (which is not our case) and, 
generally, the results are model-dependent. Therefore, in 
order to understand the influence of a complex connectiv- 
ity, we use a purely structural analysis based on the prop- 
erties of the matrix L. To this purpose, we ignore the 
intrinsic dynamics of the neurons in Eq. (0, that is, we 
just consider V = dUV . Then, there is a basis in which 
Vi ~ exp(dAit), where A; are the eigenvalues of L. It is 
well known that all the eigenvalues of the Laplacian asso- 
ciated to a network with only attractive couplings are neg- 
ative. However, when we add some repulsive connections, 
L has positive and negative eigenvalues. We find that any 
set of initial states rapidly evolves into the subspace S + as- 
sociated to the positive eigenvalues within a time smaller 
than the characteristic temporal scale of the neuronal dy- 
namics (r ~ 15 ms). 

To quantify the effect of S + , we note that, for a 
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FIG. 2: Space-time plots of the neuron voltage for a TV = 800 HH units network, with AJ = 0.2, d = 0.1, and different coupling 
connectivities: (a) Local coupling with p n = 0; (b) network with long range couplings, p — p c = 0.0055, and p n = 0; (c) p — p c — 
0.0055 and p n = 0.3; (d) p = p c = 0.0055 and p n = 0.45; (e) p = 0.015 and p n = 0.3. 
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FIG. 3: MFR (left) and network coherence cry (right) as a func- 
tion of p for several p n in a N — 800 network. Each point is 
averaged over 100 simulations, 1 s long (transients avoided), for 
different network and initial conditions realizations. Note that the 
legend applies to both figures. 



given positive A,, e dXi is a measure of how much the 
system spreads into the subspace defined by the corre- 
sponding eigenvector. Then, the ratio e dXi 1 / 'e dAmaxt = 
e d(Ai-Amax)t measures how different is the evolution in that 
subspace with respect to the one where the system de- 
velops faster. By defining the geometric average g(t) = 



n 



N p d(\ t -\ m ax)t/N 



e d{{\) A max )t ; we can es timate the 



homogeneity of the evolution in S + with a number in 
(0, 1]. A value close to 1 means the system evolves simi- 
larly in all dimensions of S + , whereas a low g implies that 
its behavior is determined by those vectors with the largest 
associated eigenvalues. 

We are interested in the behavior of g(t) as a function of 
p and p n . As the shape of g(t) with p is not very sensitive 
to time, we fix t = d^ 1 ~ r to focus our study within 
the time scale of our dynamical unit (see Fig.@lleft)). We 
observe that g = g(r) presents a minimum at p c which 
is lower for higher values of p n , and whose position shifts 
to higher p as p n increases, as observed in the numerical 



simulations (Fig. fright)). This means that, for values of 
p far from p c , i.e. where g ~ 1, the global dynamics 
is basically determined by only one positive eigenvalue, 
V(t) = V exp(Ai). On the contrary, for values of p 
close to p c , we need to consider not just one but several 
eigenvalues (the largest ones) to account for the global dy- 
namics. Therefore, the intrinsic dynamics of the system is 
minimally constrained by the structure that arises around 
p c due to the repulsive shortcuts. 

To analyze if the previous structural result affects other 
dynamics imposed on it, we consider a discrete spin-like 
dynamics in which each node i has only two possible states 
Si = dbl, with a probability p_ of being Sj = —1. Conse- 
quently, with the same link matrix L studied above, node 
i receives an input hi = LaS j € [—2, 2]. Hence, as 
other authors have pointed out 1241. 12511 . these spin-like net- 
works can be regarded as a pattern of the internal states and 
their evolution represent the global dynamics. Notice that 
the neighbor vertices linked repulsively contribute to the 
input with the opposite state. Then, in this model it is im- 
plicit that nodes linked with an attractive connection tend 
to follow the same evolution, whereas repulsive connection 
leads them to evolve differently. 

We can prove analytically that the distribution of hi 
presents two peaks: [Ax = —2p_sjl — 4j?„(l — p n ) and 

fj,2 = 2(1 — p_)wl — 4p n (l — p n ). Then, we choose to 
evolve the network according to the following local major- 
ity rule: the new state of node i is updated to S;(n + 1) = 
+1 if hi(n) > fj, 2 , Sj(n + 1) = —1 if hi(n) < \x\ and 
Sj(n + 1) = Sj(n) otherwise (i.e., the vertex keeps its 
state). 

When we average the mean state of nodes after a tran- 
sient, we find that the system changes its behavior at p « 
p c . However, since these changes in the average depend on 
p n and p_, we focus our attention in the deviation a m of 
the mean state. This quantity measures how many differ- 
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FIG. 4: Left: Dependence of g with the adding probability p, in 
a log-linear scale, for different probabilities p n . Each point is an 
average over 100 different realizations of a N = 800 network. 
Right: Deviation a m \fN of the mean state vs. p for different p n 
values in a log-linear scale, with N = 800. Each point averages 
1000 runs after a transient of 100 iterations and fixed p_ — 0.1. 

ent states are allowed for that particular network structure. 
It can be seen in Fig. lUright) that the maximum of a m is 
reached when the GCC associated to the long-range links 
spans the whole network with a minimal number of links. 

This shows how p and p n contribute to improve the syn- 
chronization even for this discrete dynamics. When the net 
is essentially a lattice, the system remains in disorder since 
there are few A^ > to spread the activity throughout. 
If we have a fully connected network with many A^ > 0, 
the whole system is activated but, since all dimensions in 
S + contribute similarly to the dynamics, the topology con- 
strains the neurons to evolve alike when they have different 
intrinsic dynamics. On the contrary, the structure close to 
p c , due to the presence of phase-repulsive links is such that, 
not only the activity is enhanced, but also the topology is 
compatible with the heterogeneity of the system. 

In summary, we have shown how a small fraction of 
phase-repulsive links can enhance synchronization in a 
complex network of dynamical units. A structural analy- 
sis allows us to obtain information about how the topology 
influences the dynamics. Surprisingly, around p c , the ver- 
satility arising from the network structure due to p n drives 
the system to a more ordered state, while far from p c the 
stiffness of the structure freezes the initial disorder. 
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